Factorial Moments in a Generalized Lattice Gas Model 
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We construct a simple multicomponent lattice gas model in one dimension in which each site 
can either be empty or occupied by at most one particle of any one of D species. Particles interact 
with a nearest neighbor interaction which depends on the species involved. This model is capable of 
reproducing the relations between factorial moments observed in high-energy scattering experiments 
for moderate values of D. The factorial moments of the negative binomial distribution can be 
obtained exactly in the limit as D becomes large, and two suitable prescriptions involving randomly 
drawn nearest neighbor interactions are given. These results indicate the need for considerable care 
in any attempt to extract information regarding possible critical phenomena from empirical factorial 
moments. 



PACS numbers: 13.85.Hd,13.65.+i,24.60.Lz,25.75.+r 



I. INTRODUCTION 

Factorial moments provide a useful tool in the analysis of high energy scattering data such as obtained in e+e" 
scattering pp scattering (at energies up to 900 GcV) |^| or the scattering of protons or heavy ions {e.g., ^^O and 
■^^S) by heavy nuclei at a projectile energy of 200 GeV/A ^. One considers the full range of some variable (usually 
the rapidity) for which one knows the average multiplicity, (n), and its dispersion, (An^) — (n^) — (n)^. The data 
is then broken into M equal bins, and one constructs the factorial moments as a suitably normalized average of a 
combination of moments in each individual bin. Specifically, 
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Although interest in these forms has been motivated by theoretical considerations Q, factorial moments are merely 
one way to arrange the experimental data. Evidently, Fq{AI) probes certain combinations of the properties of the 
1-, 2-, 3-, . . . , g-body correlations between particles in each bin. If non-interacting particles are distributed in the 
various bins in a purely statistical fashion, all factorial moments are 1 for all M. This result is completely inconsistent 
with the data. Consider, for example, pp scattering at 900 GeV. Here, F2{M) is seen to grow from 1 (for small M) 
to roughly 1.7 (for the largest values of M considered). Over the same range, F3(M) ranges from 1 to 4.2, F4{M) 
ranges from 1 to 15 and F^{M) ranges from 1 to 80. This growth in the factorial moments with M is often called 
'intermittency' in the literature and is frequently regarded as indicating the presence of fluctuations of many different 



It has been observed g that the relations between the various factorial moments for pp scattering (and all of the 
other physical processes mentioned above) can be reproduced with remarkable accuracy by the 'negative binomial 
distribution' (NB) for which 



F^^{M) = (1 + cM)(l + 2cM) ...{l + [q- l]cM) 



with 



(An2) - (n) 



(2) 



(3) 



Plots of Fq{M) versus M evidently depend on the global averages (n) and (An^). Such plots will be different for 
the various physical processes considered and are not reproduced by the negative binomial distributions as given by 
Eq. (|) §. However, since is simply (1 + cM), it is tempting to consider the more 'universal' plots of Fq versus 
F2 n|7 Such plots no longer depend on the parameter c and invite the comparison of data from very different processes 
[^. Such plots have been made over the available range 1 < F2 < 1.7. They do reveal universal behaviour |^ and are 
in striking agreement with the curves obtained from the negative binomial distribution. 
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The dramatic growth of the factorial moments with M led Bialas and Peschanski to note that 'an observation of 
a variation in (our Fq) with 5y (our M) indicates the presence of genuine fluctuations which must have some 
physical origin' Q. Many authors have taken up the challenge of describing this physical origin Some have 

concentrated on the apparent presence of fluctuations of many different sizes and considered models incorporating a 
variety of critical phenomena 1 11 1 . Others have studied both schematic and more realistic versions of cascade models 



a. 

In two interesting papers, Chau and Huang have offered a different kind of insight |13|. They imagine that the 
full range of rapidity corresponds to the N sites of a one-dimensional Ising (or lattice gas) model [|4j . This model 
is exactly solvable. They consider the g-body correlations and the factorial moments which come from a lattice gas 
model (in the limit TV ^ oo). The two parameters of the Hamiltonian are determined by fixing the global values of 
(n) and (An^). The resulting factorial moments have the form 

^rw-Y. ^^_,^l;:i^ icMf (4) 

where c is given by Eq. (|^). This result is both rather more and rather less than meets the eye. 

Let us address the 'less' first. While the lattice gas factorial moments are not identical to the corresponding moments 
of the negative binomial distribution, both share the common small M expansion 

Fq{M) = 1 + '^^'^ ~ + 0{c^M^) . (5) 

We have already noted that the leading term is due to 'one-body effects'. It should come as no surprise that the 
term of order M is precisely due to two-body correlations. Under the assumption that genuine two-body correlations 
have a finite range, the g-dependence of the term of order M is uniquely determined. The coefficient c, which can 
be expressed as a suitable integral of the two-body correlation function, is also fixed by the global dispersion, (An^). 
Thus, the observation of Chau and Huang that the factorial moments 'given by single negative binomials are almost 
exactly the same as ours for F2 1' is a trivial consequence of the rules of the game {i.e., the fixing of (n) and 
(An^)) and the fact that the Ising model predicts many-body correlations of a finite range. 

Chau and Huang further note that the negative binomial distributions are more successful than the lattice gas 
model for values of F2 larger than 1.5. This is the domain where terms of order and higher — which are not 
model independent — play a role. A few comments on the existing data are in order before continuing. Data for 200 
GeV/A ^^O and "^^S scattering on emulsions is limited to the range 1 < F2 < 1.35. Thus, this data never really probes 
the interesting (model dependent) regions of F2. (Over this range, the difference between F^'^ and F^^ is less than 
2.8% which is small compared with the experimental uncertainties.) The data for 200 GeV p on emulsion and for 900 
GeV pp scattering has a somewhat wider range covering 1 < F2 < 1.7. At the upper end of this range, the difference 
between lattice gas and negative binomial forms for F^ grows to 6.2% while the non-linear terms in the negative 
binomial form account for 24% of F3. Uncertainties in the data in this region (both systematic and statistical) are 
approximately 7.2%. The situation is similar for the existing data on F4 and F^. It is possible to find convincing 
empirical evidence for higher than linear terms in Fq. It is more difficult to distinguish empirically between negative 
binomial and lattice gas models although the former does provide a better global fit. 

We believe that Chau and Huang have provided something very positive in their work. Their model admits an 
elementary generalization which shows one possible route which leads smoothly from the lattice gas model to the 
negative binomial distributions p5[ |. While the nature of this route may be of limited empirical value given the 
current status of the data, it strikes us as being of some theoretical importance. This generalization, which is the 
primary focus of this paper, is simply stated: Consider a one dimensional lattice gas model (again in the limit where 
the number of sites, N, is infinite) where each site can either be empty or occupied by one particle which can be 
of any one of D species. Each species has a chemical potential, fi^, and each pair of species has a nearest neighbor 
interaction of strength edd' ■ Like the ordinary lattice gas (with D = I), this model allows for the analytic determination 
of the g-body correlation functions and the factorial moments, Fq{M), using the obvious generalization of standard 
techniques. As usual, this involves the construction and diagonalization of a matrix, Ai, related to the partition 
function for one pair of adjacent sites. 

The factorial moments for this generalized lattice gas model can then be specified exactly in terms of certain 
(weighted) moments of an elementary function of the eigenvalues of A4. We shall show that, for any value of D, 
these moments can be chosen to reproduce all factorial moments of the ordinary lattice gas model. Further, we shall 
show that these moments can also be chosen to approach the factorial moments of the negative binomial distribution. 
In the limit D —^ 00, this model permits the exact reproduction of all factorial moments of the negative binomial 
distribution. 
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Given the specific form of M, there is no guarantee that these conditions on its eigenvalues can be met for any 
choice of the parameters appearing in the related Hamiltonian. We shall, thus, give two methods for the selection 
of the various chemical potentials and nearest neighbor interaction parameters in A4 which explicitly reproduce the 
factorial moments of the negative binomial distribution. This choice is essentially a 'random dynamics' model in 
which the D chemical potentials are set equal and in which the various nearest neighbor interactions are drawn at 
random according to a certain distribution. The random nature of this model and the smallness of the dispersion in 
its factorial moments for D finite and small provides some understanding for the success of cascade model calculations 
of these processes and for the relative insensitivity of such calculations to the details of input parameters. 

The organization of this paper will be as follows. In Sec. ^ we summarize the results of Chau and Huang. One 
purpose of this summary is to establish the notation of the lattice gas model (as opposed to the Ising model) which 
will later be generalized. In the process we shall establish a useful (and more general) theorem regarding the factorial 
moments and shall explain precisely which features of the Ising model are responsible for the form of Eq. (4). Along 
the way, we shall point out why the general form of Eq. (5) is more a matter of definition than of physical content. 

In Sec. Ill we shall establish our generalization to a multicomponent lattice gas model containing D species, outline 
the techniques for its analytic solution, and present the general form for the resulting factorial moments. 

In Sec. 1^ we shall consider the constraints which must be imposed on this generalized model if it is (i) to reproduce 
the results of Chau and Huang or (ii) to reproduce the factorial moments of the negative binomial distribution. These 
constraints will be established analytically as a well-defined distribution of the eigenvalues of a certain {D + 1)- 
dimensional matrix. We shall also report the results of numerical studies which allow us to express these constraints 
in terms of a well-defined but random distribution of the D{D + l)/2 nearest-neighbor interaction strengths between 
the D species. 

A variety of conclusions will be drawn in Sec. 

Two short appendices are provided to deal with more technical matters. 



II. THE LATTICE GAS MODEL 



Here, we follow the arguments of Chau and Huang with the notational difference that we consider a one-dimensional 
lattice gas rather than an Ising model. We consider N microscopic sites which can have an occupation number of 
either or 1. Particles occupying adjacent sites will experience a nearest neighbor interaction of strength e. The 
system is assumed to be cyclic in the sense that particle N interacts with particle 1 as well as with particle {N — 1). 
The Hamiltonian for the system is simply 



H = e[nin2 + n2n^ + . . . + nN-iriN + riArni 
The partition function for this system is mM 



exp 



-H 



N 

1=1 



(6) 



(7) 



The external summation here covers the 2^ terms where each Ui has the value or 1. This problem is rendered trivial 
by constructing a two-dimensional matrix, M, such that 



Mmn2 = CXp 



1 



1 



-enin2 - -^ni - -^n2 



It is convenient to let the matrix indices run over and 1. 

The partition function for the system is now simply the trace of A^^. 
defining the orthogonal matrix, 0, which diagonalizes A4. 

Md^O'^Md . 



(8) 



This trace is most easily constructed by 



(9) 



The diagonal form Aid has two eigenvalues, A±, and we shall choose 9 such that {Aid)QO is the eigenvalue of larger 
magnitude, A+. (Given the form of Ai, it is clear that A+ > 0.) We then find that 



Ti[M^] = Tr[0Md"O' ] ^ Ty[M 



(10) 



In the limit 
N. 



oo, we can neglect the term A;'^ at the cost of introducing an error which is exponentially small in 
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In order to calculate the various correlation functions, we introduce the number operator matrix for site i, rii. 
Clearly, {ni)nin2 = ^in^Sin^, i-e., 




1 



(11) 



The number operator at each site is an idempotent with nf — rii. Since all sites are equivalent (due to the cyclic form 
of Eq. (H)), the average number of particles per site is simply 

= Tr[n,A^^]/Tr[A^^] . (12) 

Given the fact that Tr[A1^] — in the large N limit, we see that it is useful to replace A^d by another diagonal 
matrix, AAd, in which each diagonal element is simply divided by A+. Thus, (A^d)oo = 1 and {M.d)ii = A_ where 
|A_| < 1. With this notation, 

{n,)^TY[n,eM^e^] . (13) 

Given the trivial structure of n^, 

(n,) - (0^^0^)ii (0io)' (14) 

where the last form applies in the limit oo. While Oiq and A_ can be determined in terms of e and ^ following 

an explicit diagonalization of A^, this is not necessary. It is our intention to impose global constraints on (n) = N{ni) 
and {/S.TT?) . This can be done with equal ease working with 9iq and A_. 
The two-body correlation function can now be written immediately as 

= (0>l^X)ii(^-^r'^^)ii ■ (15) 
Taking the large N limit and making the usual assumption that N ~ j \s always 0{N), we find 

(nirii+j) = (ni)[{ni) + (1 - (ni))Ai] . (16) 

Eq. (|l^) is valid for j > 0. The structure of this two-body correlation function is instructive and very general. For 
large separations (large j), {nini+j) approaches the uncorrelated, statistical value of (rii)'^. The approach to this 
asymptotic value is exponentially fast. For small j, there are significant short-range correlations. For j = 0, we find 
{niTii^j) = (rii) independent of A_ which is merely a reflection of the fact that the number operator is idempotent. 

It is elementary to determine the global dispersion, (An^). Since it is our intention to constrain (An^) to be 
consistent with data, let us deal with this problem immediately. We first construct 

N 

= J2 ^N{n,)+2 ^ {n,,n,,) . (17) 

■il,''2 = l il<i2 

It is, of course, possible to do the summations in Eq. ( [l7| ) exactly given the form of Eq. (p^. For purposes of later 
arguments, it is more useful to obtain (n^) approximately by making approximations which neglect terms of order 
1/A^. Thus, we write 

(rO = (n)+2(n,)' l + 2(n,)(l-(n,)) ^ >^r'' . (18) 

il<i2 il<i2 

The first sum in Eq. ( p^ ) will be approximated by N'^/2. The second sum is more interesting. Since Xi_ converges 
exponentially with j, we shall allow the 12 sum to extend over all values 1 < (^2 — ii) < 00. The sum over ii then 
merely introduces a factor of N. Further, we neglect (rii) compared to 1. These approximations each introduce errors 
of order 1/A^ which are acceptable as — > cxd. We immediately obtain 

{n') = {n) + {n)^ + 2{n)^^ . (19) 

We are now able to set A_ in order to reproduce any desired dispersion. 
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The calculation of factorial moments can be simplified significantly when the number operator is idempotent. This 
issue is addressed in Appendix |^ where appropriate operator (and ensemble averag e) id entities are established. Using 
the definition of the factorial moments, Eq. ([|), using the simplifying result of Eq. ([Ad), and adopting the same large 
N approximations, we can obtain F2{M). (Now, the sums analogous to those in Eq. (|17|) extend to an upper limit of 
Nb where Nb — N/M. Since the number of bins, M, is finite, we are also concerned with the limit Nb oo.) We 
find 



F^^{M) = l + cM 



(20) 



where 



(n) 1 - A_ 



Returning to Eqs. (|l|) and (^, it is now elementary to construct the g-body correlation functions as 



I'll ' +22 



= + (1 - {n,))Xl][{n,) + (1 - (n,))AL^] . . . [(n,) + (1 - {n,))\': 



(21) 



(22) 



Eq. (22) applies for ^l,^2^ ■ ■ ■ ,iq > 0. This form, which is not unique to the lattice gas model, determines the form 
of the factorial moments and merits comment. The g-body correlation function is completely determined by the 
two-body correlation function and (rii). More precisely, the g-body correlation function is a product of the {q — 1) 
two-body correlators between adjacent particles. Finally, the individual two-body correlators are given as the sum 
of a statistical term and a short-range piece. These qualitative features are sufficient to determine the associated 
factorial moments uniquely without any detailed information about the precise nature of the short-range part of the 
two-body correlator. To emphasize this independence, we shall make the replacement of (1 — {ni))Xi_ by g{j) where 
the only restriction on g{j) is that it is of short range. 

Now let us turn to an arbitrary factorial moment in the limit N ^ oo. We again use the operator identity derived 
in Appendix |A|, Eq. ( |A6| ), and replace sums by integrals to obtain 



LG 



(M) = q\ 



M 



Nb 



dxi 



Nb~xi 



dX2 



Nb-xi- 



dXn 



[7^+5^)] 



(23) 



The structure of this result is clear. We now imagine expanding the {q — 1) square brackets. Since the function g{xi) 
is short range, we can extend the Xi integration to infinity for any of those 2'^~^ terms which contains a factor of 
g{xi). Let us introduce the notation 



G = 



dx g{x) 



(24) 



The term 0(Af°) can only come from picking the statistical factor in each term. There is 1 way to do this. We obtain 
no factors of G. The remaining integrals give a factor of Ng/ql. Hence, the leading term in F^'~^{M) is always one, 
as expected. The term 0{M) comes from picking (g — 1) statistical factors and 1 factor of g. There are {q — 1) ways 
to do this, and we obtain a factor of G. The remaining integrals give a factor of iV^^ /(<? ^ !)!■ Hence, the next term 

l)/2. The general result is now obvious: 



in F^^iM) is precisely {2GM/{n))q{q 



9-1 r 



F^^^iM) = y: 



k=0 



2GM 



q\{q-iy. 



{q-k)\k\{q-l-k)\2^ ■ 



Comparing Eq. (25) with q — 2 with the result of Eq. (pO[), we see that we can make the identification 

2G 



(25) 



(26) 



Eqs. ( p5| ) and (|2^) are precisely the forms found by Chau and Huang as quoted in Eq. (|4|). The present derivation 
makes it clear that these results do depend on the fact that the g-body correlator is a product of the (q — 1) consecutive 
two-body correlators and that they do not depend on the specific form of the short-range two-body correlations. 

The fact that 
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Fq{M) = 1 



cM + 0{c^M'^) 



(27) 



is even more general. It depends only on the fact that the g-body correlation function contains no long-range terms 
and thus includes all one-dimensional models. For any model we can consider 



+i2 



'■11+. ..+1, 



(28) 



in the limit where one of the indices {e.g., ik) is small and all others are large. In the absence of long-range correlations, 
this limit must reduce to 



(29) 



(Translational invariance ensures that the final term here depends only on ik-) Such terms are the only ones which 
contribute to the linear piece oi Fq{M). Their q dependence is protected by elementary combinatorics. The coefficient 
c is protected by the fact that the two-body correlation has been adjusted to fit the empirical global dispersion, (An^). 
In a similar fashion, the terms O(Af^) receive contributions from both two- and three-body correlations. If we had 
insisted on establishing global constraints on (n^), these terms would be similarly protected. 

Thus, the fact that lattice gas factorial moments agree with those coming from the negative binomial distribution 
in the limit cM —^ provides no indication that the lattice gas model is correct. It is merely an indication of the 
fact that there are no long-range correlations and that the 'rules of the game' are to consider only models with fixed 



global 



and (At 



III. A GENERALIZED LATTICE GAS MODEL 



As noted in Sec. ^ the ordinary lattice gas model does not reproduce the factorial moments of the negative binomial 
distribution. The fact that there is agreement in the small cM limit is a consequence of general properties of one- 
dimensional models and is not indicative of any particular merit of the lattice gas model. In this section we shall 
construct a simple generalization of this lattice gas model which has the capacity to reproduce the negative binomial 
factorial moments exactly. This model is a conceptually simple one-dimensional model which is also exactly solvable 
(in a sense which is appropriate for constructing factorial moments). 

Consider a lattice gas in which site i is either unoccupied or is occupied by at most one particle which can be of any 
one of D species. Each species has its own chemical potential, fid- Particles again have a nearest neighbor interaction. 
The strength of the interaction between a particle of species d and a particle of species d' is e^d' ■ 

This model can be solved using precisely the standard techniques of Sec. The matrix A4 now becomes a real, 
symmetric matrix of dimension (D + 1). The elements in this matrix are 



A^oo = 1 

Mm = cxp [-Hd/2] 
Mdd' = cxp [-edd' - A*<i/2 - Hd'/'i] 



(30) 



(We allow the matrix indices to run from to D. Again, the inverse temperature, f3, has been set equal to 1.) The 
number operator at site i is also a {D + l)-dimensional matrix having the form rii = 1 — T with Tdd' — SdoSd'o, i.e., 







V 



(31) 



1/ 



Both this number operator and T are idempotent, so that the results of Eqs. (A3) and (A6) remain valid both at the 
operator level and at the level of ensemble averages. 

We again find the transformation 9 which diagonalizes M and places the eigenvalue of largest magnitude in (A^d)oo 
[0. As before, we construct the matrix by dividing each element of A^^; by the eigenvalue of largest magnitude. 
(The form of again ensures that this largest eigenvalue will be positive.) Thus, {M.d)ovi is again 1. While we shall 
retain this notation, it is not completely necessary in practice. The spirit of the model is, ultimately, to take the limit 
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as A'' ^ (X) for fixed (n). In this hmit, the average occupation per site, (rii), tends to zero. This limit will be realized 
by taking the limit as all of the /z^ +00. In this limit, the coupling of the D x D submatrix (with d, d! 7^ 0) to the 
remaining elements of M. becomes arbitrarily weak. In the limit, it is legitimate to treat this coupling in lowest-order 
perturbation theory. We shall return to this point in Sec. 



IV. 



With these preliminaries in hand, the construction of the various correlation functions and factorial moments 
proceeds as before. The average number of particles per site is simply 

{n,)^^A = Tr[{t-T)eM^9^] (32) 

In the limit N ^ 00, Ai^ may be set equal to T with exponentially small errors. The two terms in the factor (1 — T) 
must be treated separately [Q. One immediately finds 

(n,> - 1 - {9nof (33) 

As usual, (n) and hence (^oo)^ will be fixed by experiment. 
The two-body correlator is readily calculated as 

= Tr[{l - T)0M'^9'^il - T)9T9^] . (34) 

Expanding the factors of (1 — T), one immediately obtains 

D 

(nin.+j) = (n,)[(n,) + (l-(n,»^((?o^)'P,] . (35) 

£=1 

This form is virtually identical to that of Eq. (p^. The first term is again the asymptotic statistical probability of 
finding one particle (of any species) on each of the two sites. Since |A^| < 1, the second term here again vanishes 
exponentially for large j. However, in this case, we are confronted with a sum of exponentials rather than a single 
term. 

Following the arguments from Eq. (|l^) to (|2l| ) we find that ||l9| 



F^'^'^ (M) = I + cM (36) 



with 



_ {n^)-{nf-{n) 2 ^ A, 
{ny (n) ^ 1 - 

Here, we have taken the notational liberty of exploiting the fact that 9 is an orthogonal matrix and made the 
substitution 

iOoi? ^.^^ j:bj = l . (38) 

Again, the similarities between Eqs. (|3^) and (|l]) are striking. The only difference is the presence of the sum over £. 

Since certain differences arise at the level of the three-body correlation function, this term is worth discussing 
specifically. Obviously, 

(n,n,+,n,+j+fe) = Tr[{l - T)9M'Xi'^ " T)9M^Xi'^ " T)9T9^] (39) 
Expanding the factors of (1 — T) and constructing the traces gives rise to the form 

+ (n,)(l-(n,))^6|Af'= (40) 

This result differs from the expression for the ordinary lattice gas. The three-body correlator cannot be expressed 
as a product of consecutive two-body correlators. This is important since it will give us precisely the freedom we 
need to proceed from the lattice gas factorial moments to those of the negative binomial distribution. Nevertheless, 
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this correlator is still determined by the two-body correlation function. It merely requires the addition of (nirii^j^k)- 
This fact is also important, since we seek factorial moments which depend on the global properties of the distribution 
through c only. Note that in the special case of D = 1 appropriate for the ordinary lattice gas, the sums collapse to 
a single term with — 1. 

It is now straightforward to perform the sums necessary to determine F^'"^{M) in the limit as — > (X). Using 
Eq. (^) we find 



^^^{M) = l + 3cM + -(fM^ . 



(41) 



The presence of the term involving {nirii-f-j^k) in the three-body correlation function has forced the introduction of 
the new term 



D 



A. 



(42) 



We note that this third factorial moment will equal that of the original lattice gas model provided that cP = . There 
are (at least) two situations where this will happen. One is the case where precisely one of the 6^ is equal to 1 (with 
the remaining 6^ equal to 0). The other is the case where all of the \i are equal. 

We now turn to the general case of F^^'^(M). The strategy is now perfectly straightforward and exceptionally 
tedious. We construct the g-body correlation function as the obvious generali zati on of Eqs. ( ^ ) and (|3^). We expand 
the various factors of (1 — T) and perform the requisite traces. Using Eq. (A6), we then perform the requisite bin 
sums. The final result follows from some algebra and combinatorics. We find 



9-1 



fc=0 



where 



Jk — 



k\{q-k)\ 



1(9- fc) 



(43) 



(44) 



z=0 



Here, we have introduced the definition 



' 1=1 ^ 



(45) 



which implies si = c/2. 

It is also desirable to derive a generating function from which the factorial moments can be obtained by simple 
differentiation. To this end, we first observe that flf' can be re-written as 



Aq) _ 

Jk — 



1 



{q-k)\ 



dzi ^ ^ 



(46) 



-I z=0 



We now change the summation index in Eq. ( ^ ) from fc to (g — fc) to obtain 



^1 



d1 

dzi 



9 / oo \ 

k=l \ ti=l / 



(47) 



2 = 



The sum over k can be extended to range from to c» since the corresponding derivatives are zero. The factorial 
moments thus become 



f: 



GLG 



di 



(48) 



z=0 



Eqs. 



-(Eq) and (Es) represent the final results of our generalized lattice gas model. 
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The various terms, fl , which characterize the factorial moments of our generalized lattice gas model are completely 

determined by the leading terms for k' < k. This follows from combinatoric arguments and has virtually 

nothing to do with the underlying microscopic details of the model. It does involve two basic ingredients of the 
model. First, that each site has an occupancy of either or 1. Second, that the non-statistical pieces of the various 
correlations are of short range [ pO| . It was not a priori obvious that any one-dimensional model would yield the 
factorial moments of the negative binomial distribution. It is now clear, however, that the generalized lattice gas 
model is capable of reproducing the factorial moments of any distribution (including the NB) through a suitable 
choice of the s^. It remains to be demonstrated that there exists a choice of the underlying Hamiltonian which will 
provide the desired Sfj,. Such demonstration is the topic of the next section. 



IV. OBTAINING THE FACTORIAL MOMENTS OF THE NEGATIVE BINOMIAL DISTRIBUTION 

As we have repeatedly emphasized, we are certainly able to reproduce the factorial moments of the ordinary lattice 
jas model. We will show in Appendix ^ that this happens in the special case where 

As noted, this special case can be realized either when the sums over i contain only one term or when all the reduced 
eigenvalues are equal. 

We are also able to pick the various terms, s^, in such a way as to reduce Eqs. (p^-(p5|) to Eq. (||) exactly. It is 
readily verified that the choice 

(50) 



/i+ 1 



establishes the desired equality. (The demonstration of this fact is given in Appendix ^.) With this choice, the 
factorial moments of our generalized lattice gas model are rendered identical to the factorial moments of the negative 
binomial distribution. This answer poses a question. Is it possible to pick the dimension D (the number of species) 
and to find a nearest neighbor Hamiltonian and its related matrix, Al, such that the constraints of Eq. (^0|) are 
satisfied? As we shall see, it is. 

As noted in the previous section, we are ultimately interested in the N 00 limit. We shall arbitrarily set all 
chemical potentials equal to fio and realize this limit by taking fj,a 00. It is now sufficient to consider the D- 
dimensional submatrix, A^, obtained by neglecting the 0-th row and column of A4. Using first order perturbation 
theory, the coupling of M to the remaining elements of Ai can then be treated exactly. This will result in a slightly 
different expression for s^. We obtain 



D 



A*. 



. / I ~ 



(51) 



where the Xk are the eigenvalues of Ai |21|] and the are normalized coefficients which follow from the eigenvectors 
of M as 

1 ^ 

ak^N- -J2e,k ■ (52) 



2 



Here, is the matrix that diagonalizes M and N is chosen such that ^ 
We now re-write Eq. (BOf) as 




(53) 



M+ 1 



The purpose of this re-writing is to emphasize that, in order to reproduce the negative binomial factorial moments, 
it is necessary to satisfy an infinite number of constraints, Eq. (|5^). There are only a finite number of parameters 
in the Hamiltonian for finite D. Specifically, there are £)(£> + l)/2 interaction parameters to be specified 12^. The 
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satisfaction of an infinite number of constraints will evidently require the limit D — s- cx). (We shall return to this point 
below.) 

It is interesting to note that in the Z? — s- oo limit, we can replace the summation in Eq. (^3|) by an integral in order 
to obtain 



+1 



dX a{Xy 



A 



1 - A 



(54) 



where the coefficients a(A) now play the role of a (non-negative) continuous weighting function. Working with a new 
variable 



A 



V = 



1-A 



we thus seek another weighting function, W{r]), such that 

f + CXD 



1/2 M+1 



(55) 



(56) 



As it happens, this set of constraints allows us to determine W{r]) by inspection. We simply obtain a step function, 

W{tj) ^ 0{c{n) - r])0{7j) . (57) 

This choice (barring questions of uniform convergence) is unique. It is readily transformed into a statement about 
the weighting function a{Xf' . We find a(A)^ — except 



a{\f 



1 



c(n)(l- A)2 



for < A < 



c{n) 
1 + c{n) 



(58) 



This statement of the constraints will be extremely useful in finding a scheme for picking the parameters in the 
Hamiltonian in order to reproduce the negative binomial moments. 

It is clear that the form of M as given by Eq. (30) (for c?, d' ^ 0) is somewhat special. It is necessary to demonstrate 
that there exists at least one scheme for picking the D{D + l)/2 parameters of M in such a way that the various 
constraints for obtaining the negative binomial distribution are satisfied. Here, we shall simply propose two different 
prescriptions and demonstrate by numerical example that these prescriptions indeed yield the desired result. Of 
course, we shall make no claims regarding the uniqueness of our prescriptions. The first prescription is as follows: 

Draw D random numbers, Xd, from the interval [0,1]. Set M.dd equal to XdL where L > 0. Set all 
off-diagonal elements of M equal to zero. For each draw, choose L such that si = c/2. (This sets the 
dispersion to its empirical value for each draw.) 

The physical content of this prescription is clear. While identical particles experience a nearest neighbor interaction 
which ranges from « — /Xq to +00, inequivalent particles experience a nearest neighbor interaction in the range 
—fjLo <SC Cdd' ^ +00. The success of this prescription is guaranteed by a comparison of Eqs. (|5^) and (|5^). Choosing 
Jvl to be diagonal guarantees that the various sums in Eq. (|5^) will be independent of k. The A-dependence of 
Eq. ( ^8| ) will be obtained from Eq. (|5^ ) provided that the eigenvalues are distributed uniformly over a finite, positive 
range. This is ensured by our prescription. 

For our numerical studies, we consider the values of (n) — 20 and (An^) = 110 appropriate for the description 
of pp scattering at 200 GeV. (For this reaction, empirical values for the factorial moments are readily available. 
Qualitatively similar results are obtained for other high-energy scattering experiments.) These data indicate that we 
should set 



c(n) = 



(An2 



4.5 



(59) 



We now construct the ratios 



fc=i "fc 



(li)' 



(60) 
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Our only constraint, Si = c/2, ensures ri = 1. When we have satisfied the constraints of Eq. (|5^, all of these 
ratios should be equal to one |23|| . Since there are no parameters to adjust, this is an extremely stringent test of 
our prescription. It succeeds. In Table | we report the results of our numerical studies for 2 < g < 6 and D = 1, 
2, 4, ... , 512. The notation '((r^))' is a reminder that, since our theory is now randomly drawn, we must perform 
an 'ensemble average' over theories as well as the usual thermodynamic ensemble average. For each value of D, this 
ensemble average was performed over 10^ theories. 

Several comments are in order. Our prescription meets the remaining conditions of Eq. (|5^) with increasing accuracy 
as D ^ oo. For fixed q, the value of {{vq)) approaches 1 like 1/D a.s D ^ oo. The dispersion also vanishes (like 
1/a/D). Thus, as D becomes large, our simple prescription converges to the results of the NB for any fixed q. For 
fixed the error in {{rq)) and its dispersion grow as g — > oo. Our point here is that there exists at least one simple 
prescription for satisfying Eq. (^0|) exactly. 

Given the limited scope and quality of existing data, we do not need particularly large D to fit the relations between 
measured factorial moments. The value of I? = 16 results in sufficiently small errors and dispersions that randomly 
drawn dynamics have a high probability of reproducing the relations between empirical factorial moments for pp 
scattering at 200 GeV within existing experimental uncertainties. This value of D = 16 is also sufficient to provide 
a quantitative description of relations between factorial moments for 900 GeV pp scattering and, indeed, of all other 
high-energy scattering experiments for which factorial moments are known. To demonstrate this. Fig. |^ shows the 
factorial moments F3 to F5 for pp scattering at 200, 546 and 900 GeV as a function of ^2- We have suppressed the 
model-independent constant and linear pieces of these moments, i.e., we plot 

F,(F2)=F,(F2)-l-fc^(F2-l) . (61) 

The theoretical 'bands' shown on each figure correspond to the theoretical dispersion obtained with D = 16. (More 
than 50% of all randomly drawn theories lie inside this band.) These bands were obtained with parameters chosen 
to reproduce the 200 GeV data. The use of parameters chosen to reproduce the higher energy data would result in 
somewhat narrower bands that are slightly shifted downwards. 

The error bars in Fig. |^ represent the published uncertainties and have roughly equal statistical and systematic 
components. Some comment regarding the statistical error is in order. Since F2{M) and Fq{M) are drawn from the 
same data set, significant statistical correlations can exist. Unfortunately, it is impossible to make a quantitative 
statement about such correlations without a detailed investigation of the data reduction techniques actually used. We 
have, however, performed simulations based on a suitably modified form of the negative binomial distribution. The 
results suggest that the uncertainties shown in Fig. |^ actually represent a significant overestimate of the real uncer- 
tainties in the determination of Fq[F2). This belief is also supported by the observation that the data points appear 
to have a far smaller spread than their errors would indicate. A reduction of these errors by a factor of approximately 
3 would make it possible to provide a convincing discrimination between the negative binomial distribution and the 
Ising model distribution. Our simulations suggest that an error reduction of this magnitude is not out of the question. 

This first prescription has the virtues of simplicity and guaranteed success. It can, however, be criticized on the 
grounds that it has singled out the species-diagonal interactions, edd, for special treatment. We thus consider a second 
and more democratic prescription in which all pairs of species are treated equivalently: 

Draw the various terms edd' at random over the uniform interval 

Co - A^o < f-dd' < Co - ^^o + Ae . 



For each choice of 7W, adjust in order to reproduce the desired value of c. (This can always be done.) 
Adjust Ae to reproduce the value r2 = 1. In the event that this cannot be done, discard the draw [p^ . 

This more democratic selection of Ai means that we must actually perform the requisite diagonalizations. This means 
that one must be content with studying a smaller number of samples. More seriously, it is no longer elementary to 
establish rules which guarantee the desired equivalence between Eqs. (^2|) and (^8|). It is for this reason that we have 
imposed the somewhat artificial constraint that r2 = 1. 

The results of numerical studies with this second prescription are summarized in Table ^ for D — 32, 64, 128 and 
256. For each value of D we considered 10^ matrices. In this case, we plot the average values (and the corresponding 
dispersions) for Co and Ae as well as the various ratios rq for 3 < g < 6. The large values of Ae required indicate 
that we are dealing with sparse matrices. This has two apparently negative consequences. First, there will be a great 



many eigenvalues close to zero. Second, there will be only a small preference for positive eigenvalues |25|. These facts 
would appear to make it difficult to satisfy Eq. (|5^). At the very least, we anticipate much slower convergence (as 
a function of D) when using this second prescription. On the other hand, the eigenvectors associated with negative 
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eigenvalues necessarily involve terms of mixed signs so that there will be a greater tendency for cancellation in the 
sums over 0ik in Eq. (^2|) for negative eigenvalues. 

Somewhat surprisingly, these effects conspire to a remarkable degree. The numerical results of Table || suggest 
that our second prescription also leads to the factorial moments of the negative binomial distribution in the large D 
limit. This claim of success is necessarily somewhat tentative. The slower convergence of this prescription suggests 
the need to study values of D larger than those in Table |. The need to diagonalize many large matrices restricts us 
to smaller values of D and smaller statistical samples. We have performed other tests which also suggest that this 
second prescription works. For example, we have made histograms of the of. of Eq. as a function of the variable 
ry defined above. Except for the expected 'Gibbs phenomenon', we find remarkable agreement with the distribution 
of Eq. ( ^7| ) as D becomes large. We have also studied the cases D = 512 and 1028 with very poor statistics. These 
results are also consistent with our assertion that the factorial moments of the negative binomial distribution will 
result in the large D limit. 

The purpose of this second prescription has been to demonstrate the existence of at least one democratic way to 
reproduce the factorial moments of the negative binomial distribution. More efficient schemes may well exist. 

V. DISCUSSIONS AND CONCLUSIONS 

We have demonstrated that a simple extension of the one-dimensional lattice gas model to include nearest neighbor 
interactions between pairs of D species of particles can reproduce the factorial moments of the negative binomial dis- 
tribution exactly provided that a number of constraints are met. We have further demonstrated that these constraints 
can be met by simple models of the microscopic dynamics in which all particles have the same chemical potential and 
in which the nearest neighbor interactions are drawn at random according to the prescriptions of the previous section. 

Since our model approaches the negative binomial distribution as D — s- cx), it will reproduce the relations between 
empirical factorial moments including those obtained in pp scattering, e+e^ scattering and relativistic heavy ion 
collisions. The only question is whether the number of species required to provide an adequate fit is sufficiently small 
to be considered 'physically reasonable'. Given the results of Fig. |l| for 16 species, we believe that the present model 
is physically reasonable. These results also offer some understanding for both the success of cascade calculations and 
the anecdotal observation that the results of calculations are often surprisingly insensitive to the details of the model. 
According to our picture, any randomly drawn theory would be likely to do as well (at least at the level of the factorial 
moments). 

It has been suggested that the slopes, aq/a2, in plots of InFg versus lnF2 can be used to distinguish between 
cascade systems and systems exhibiting some critical phenomenon such as the transition to a quark-gluon plasma . 
It has been shown that systems at the critical temperature of a second-order phase transition display monofractal 
structure with aq/a2 ^ {<l~ 1) PT|] . In contrast, cascade models lead to a multifractal structure with cnqja^ ~ 9(9— 1) 
(in Gaussian approximation) |12| ]. Our model (along with the Ising model) shows that the monofractal structure of 
a critical phenomenon can also be obtained from a simple but heterogeneous (equilibrium) system. It indicates that 
some caution is required before claiming insight into the nature of microscopic mechanisms on such a basis. Indeed, 
the present model may be a useful test case for proposed schemes for the detection of critical phenomena from the 
study of, e.f?., factorial moments. 

Before closing, we would like to offer one speculative remark regarding the existence (or possible non-existence) of 
signatures which would honestly discriminate between 'interesting' critical systems and those 'uninteresting' hetero- 
geneous systems considered here with which they can be confused. We suspect that the study of global properties 
— such as the factorial moments — can only provide convincing tests of criticality for systems whose time evolution 
can be studied. The present results suggest that it may always be possible to devise equilibrium, heterogeneous 
models which can simulate the manifestations of critical phenomena to arbitrary accuracy in systems where the only 
information available is that provided at t = 00. Evidently, this class of systems includes all scattering experiments 
in nuclear and high-energy physics. The only way to distinguish between these classes of models is to provide inde- 
pendent arguments that the heterogeneous alternatives to critical systems are 'too artificial' or 'too unphysical' to be 
accepted. Given the present results, we suspect that this case cannot be made. Thus, in physics as in the rest of life, 
much more information is to be gained from the study of 'living' systems than from performing autopsies on 'dead' 
systems. 
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APPENDIX A: SIMPLIFYING FACTORIAL MOMENTS 



Consider any model for which the number operator at any site i is an idempotent: 

2 



(Al) 



Such models include either the ordinary lattice gas model of Sec. |l| or the generalized lattice gas model of Sec. Ill 
Define the number operator for a single bin, 



where i runs over every site in the bin. 

Now assume that we know that the following operator identity holds for some fixed order q: 



n{n — l){n — 2) . . . (n — g + 1) = 



"^12 ^ia ■ ■ ■ "^iq 



(A2) 



(A3) 



The sum exten ds o ver all sites in a given bin. The prime denotes that no two of the site indices, ii to iq, are equal. 
Note that Eq. ( [A^ ) applies at the operator level before ensemble averages are carried ou t. O ur aim is to show that, if 
Eq. ( [A3| ) applies for q, it necessarily applies for q + 1 a,s well. To see this, multiply Eq. (A3) by the number operator, 
n. This leads to 



i[n{n - l){n - 2) . . . (n - g + 1)] = ^ 'n.^i 



(A4) 



The second term here simply corresponds to the fact that the index iq-^-i can be equal to any one of the q indices ii 
to iq. Using Eq. (A3), we can move the final term in Eq. (|A4|) to the left of the equation to obtain 



n{n - l)(n - 2) . . . {n - q + l){n - q) = ^'"m 



(A5) 



This is the desired extension of the result. Thus, Eq. (A2) will apply for all q if it applies for q = 1. It applies triv ially 
for q = I since the primed sum is of no meaning in this case (there being only one term in the product) and Eq. (A3) 
collapses to the definition, Eq. (A2). 

Given the form in which the various correlation functions can be expressed in our models, one final rearrangement 
of Eq. (A3) is useful. Specifically, 



n(n - l)(n - 2) . . . (n - g+ 1) = g! ^ 



(A6) 



il<i2<---<iq 



The final result is that Eqs. ( |A3D and (A6) apply for all q for all problems in which the number operator is 
idempotent. Further, one can take ensemble averages of these operator identities, and the ensemble averages can 
be taken under the summations. These results provide a considerable technical simplification in the evaluation of 
factorial moments. 



APPENDIX B: SPECIAL CASES OF THE GENERALIZED LATTICE GAS MODEL 



Here, we wish to reduce the factorial moments of our generalized lattice gas model to (i) the factorial moments of 
the ordinary lattice gas model and (ii) the factorial moments of the negative binomial distribution by simplifying the 
general expressions, Eqs. (^)-(^, with the use of special assumptions regarding the choice of the terms s^. 

Consider, first, the ordinary lattice gas results. The purpose here is an exercise in method since we already know 
the result to be true. Set the various as 
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We can then perform the sum in the square brackets in Eq. ( ^ ) exactly to yield 



zc/2 



(Bl) 



(B2) 



The desired derivative is immediately obtained as 



dz'' VI - zc/2 
This immediately gives 



(q-k) 



(q - - 1)! 



Aq) 
Jk 



(g-fc)!/j!(g-l-/c)!2fc 



(B3) 



(B4) 



Comparison with Eq. ^ reveals the expected agreement. 

For the case of the negative binomials, it is more convenient to start with the generalized generating function of 
Eq. (^). We now wish to set 

(B5) 



M+ 1 



The sum over /i can again be performed exactly to yield 

„ In [1 — zc\ 
1 + V z^'s^ = ^ ^ 



(B6) 



Substitution of this expression into Eq. (^) yields 



dzi 
(cM)« 



exp 



ln[l 



cM 



z=0 



d"^ 
dyi 



v=o 



{cMf 



cM J \cM 



cM 



= (1 + cM)(l + 2cM) . . . (1 + [q - l]cM) 
which is seen to be identical with the NB result of Eq. (0) . 



(B7) 
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FIG. 1. Plots of Fq as defined in Eq. (|6l|) versus F2 for 3 < g < 5. Data are taken from Ref. [g]. The dashed line corresponds 
to the predictions of the ordinary lattice gas (Ising) model. The bold face line represents the negative binomial distribution. 
The shaded area is the prediction of the generalized lattice gas model according to our first prescription for D = 16. This 
'band' around the NB becomes narrower as D increases. The error bars shown are representative for the quoted uncertainties 
in the data (see the text for further comments on the error bars). 



TABLE I. The ensemble average and dispersion of rq according to our first prescription for 2 < g < 6 and various values of 
the number of particle species, D. The average was taken over 10^ randomly drawn theories with c(n) — 4.5 for each draw. 
The case D = 1 corresponds to the ordinary lattice gas (Ising) model and has no dispersion. 



D 




((^3» 


(^4» 


(^5» 


((^6)) 


1 


0.75 


0.5 


0.3125 


0.1875 


0.109375 


2 


0.812 ±0.023 


0.598 ± 0.035 


0.417 ±0.037 


0.280 ± 0.033 


0.184 ± 0.027 


4 


0.881 ± 0.046 


0.723 ± 0.078 


0.569 ±0.091 


0.435 ± 0.091 


0.326 ± 0.083 


8 


0.941 ± 0.067 


0.851 ±0.132 


0.751 ±0.179 


0.652 ±0.207 


0.558 ±0.218 


16 


0.981 ±0.076 


0.951 ±0.173 


0.916 ±0.270 


0.879 ± 0.360 


0.842 ± 0.441 


32 


0.998 ± 0.069 


0.999 ±0.170 


1.005 ±0.295 


1.017 ±0.442 


1.039 ± 0.612 


64 


1.002 ±0.050 


1.007 ±0.128 


1.019 ±0.229 


1.040 ± 0.355 


1.071 ±0.515 


128 


1.001 ±0.035 


1.005 ±0.088 


1.013 ±0.154 


1.026 ±0.233 


1.044 ± 0.327 


256 


1.001 ±0.024 


1.003 ±0.061 


1.007 ±0.105 


1.013 ±0.155 


1.023 ±0.210 


512 


1.000 ±0.017 


1.002 ± 0.042 


1.004 ± 0.072 


1.007 ±0.106 


1.012 ±0.142 



TABLE II. The ensemble average and dispersion of Vq according to our second prescription for 3 < g < 6 and various values 
of the number of particle species, D. The average was taken over 10^ randomly drawn theories with c(n) = 4.5 for each draw. 



Also shown 


are average 


value and dispersion 


of to and Se. 








D 


to 


Se 


((^3)) 


(('-4)> 


{{rs}) 


((^6)) 


32 


0.21 ±0.37 


114 ±75 


0.945 ± 0.025 


0.862 ± 0.059 


0.769 ± 0.094 


0.674 ± 0.124 


64 


0.47 ±0.23 


149 ± 63 


0.959 ±0.022 


0.895 ± 0.054 


0.820 ± 0.089 


0.741 ± 0.122 


128 


0.62 ±0.18 


223 ± 71 


0.973 ± 0.021 


0.930 ± 0.054 


0.876 ± 0.093 


0.818 ±0.133 


256 


0.73 ±0.15 


363 ± 85 


0.988 ± 0.020 


0.967 ± 0.054 


0.941 ±0.098 


0.910 ±0.148 
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